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NOMENCLATURE 


Symbol 


9 

k 

P 

Po 

AP 

Q 

R 

T 

AT 

u 

v 

X 

AX 

y 

Ay 


specific  heat,  [J  kg"1  K_1] 
arbitrary  value,  eq . (15 ) 
acceleration  due  to  gravity,  [ms-1] 
thermal  conductivity.,  [W  nr1  K_1] 
pressure,  [Pa] 

pressure  at  original  point,  = R % T0 

pressure  difference,  = P - P0 

total  heat  is  added  into  the  enclosure,  [W] 

gas  constant,  [J  kg-1  K"1] 

temperature,  [K] 

temperature  difference,  = Tw  - T0 
velocity  component  in  x direction,  [m  s"1] 
velocity  component  in  y direction  [m  s_1] 
coordinate  along  the  vertical  wall,  [m] 
grid  mesh  size  in  the  x direction,  [m] 
coordinate  along  the  horizontal  wall,  [m] 
grid  mesh  size  in  the  y direction,  [m] 


Greek  Letters 


& 

y 

P 

p 

r 

At 


coefficient  of  volumetric  expansion,  [K_1] 

specific  heat  ratio,  Cp/Cy 

absolute  viscosity,  [Pa  s] 

density,  [kg  m~3] 

time,  [s] 

time  step,  [s] 


Subscripts  and  Superscript 


0 

1 
j 

m 

n 

r 

w 


value  at  the  wall  except  the  heated  wall 
subscript  denoting  the  i th  grid  point  in 
subscript  denoting  the  j th  grid  point  in 
average  value 

superscript  denoting  the  time  at  Tn 

restricted  value 

value  at  the  heating  surface 


the 

the 


i i i 


di recti  on 
di recti on 


COMPUTATION  OF  TWO-DIMENSIONAL  TIME-DEPENDENT  NATURAL  CONVECTION 
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National  Engineering  Laboratory 
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Studies  of  natural  convection  processes  generally  assume 
an  incompressible  fluid  wherein  the  density  is  a function  of 
temperature  only  (the  Boussinesq  approximation).  However, 
local  pressure  gradients  caused  by  rapid  variations  in  the 
heated  wall  temperature  cannot  be  described  within  this  ap- 
proximation. These  time-varying  gradients  cause  fluid  motions 
which  perturb  the  quasi-static  natural  convection  process.  In 
this  study,  we  describe  a numerical  analysis  procedure  which 
includes  compressibility  effects  and  allows  computation  of 
transient  fluid  motions  during  onset  of  natural  convection. 
Details  of  the  computational  procedure  and  preliminary  results 
for  one  geometry  are  given. 

Key  words:  compressible  fluid  motion;  convection;  finite 
difference  approximation;  heat  transfer;  natural  convection; 
nonlinear  convection;  numerical  integration;  transient  fluid 
motion;  transient  heat  transfer. 
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INTRODUCTION 


For  two-dimensional  time-dependent  laminar  natural  convection  about  a 

/ p O C Cl 

heated  surface,  many  numerical  solutions  are  known  in  the  literature.  5 ’’ 
These  solutions  use  Boussinesq's  approximation  in  which  the  thermophysical 
properties  are  taken  to  be  constant  and  the  fluid  is  assumed  incompressible, 
except  when  considering  the  body  force  term  in  the  equation  of  motion.  These 
assumptions  lead  to  a numerical  formulation  in  which  the  pressure  terms  are 
eliminated  from  the  equation  of  motion  and  a stream  function  is  defined  which 
satisfies  the  equation  of  continuity.  Spiegel  and  Veroni's  (4)  have  presented 
the  conditions  under  which  the  Boussinesq  approximation  is  applicable  for 
thermal  convection  in  compressible  fluids;  although  the  Boussinesq  approximation 
is  valid  for  a number  of  natural  convecting  problems,  it  is  of  uncertain 
accuracy  for  studies  of  supercritical  fluid  motion  where  the  thermophysical 
properties  of  these  fluids  are  strongly  dependent  on  temperature  and  pressure. 
However,  pressure  changes  generated  by  pulsed  thermal  input,  which  are 
significant  in  determining  the  fluid  motion  in  early  time  periods,  cannot  be 
described  within  such  a formulation. 

From  this  point  of  view,  the  authors  have  attempted  to  get  computational 
results  using  the  program  PDETWQ  (7)  for  two-dimensional  time-dependent  natural 
convection  of  a compressible  fluid  in  a rectangular  enclosure.  However,  stabil- 
ity problems  arose  in  computational  work,  perhaps  associated  with  implicit  time 
integration.  Also,  the  cross  differential  term  d^/{dxdy)  in  momentum 
equations  is  neglected  in  the  PEDETWO  program,  and  it  was  not  known  whether  this 
could  lead  to  a significant  error. 

In  this  paper,  computational  analysis  is  described  on  the  laminar  natural 
convection  heat  transfer  from  an  isothermal  wall  to  compressible  fluid  within  a 
rectangular  enclosure,  taking  into  account  the  variation  of  thermophysical  pro- 
perties and  cross  differential  terms  in  momentum  equations.  Calculated  results 
for  both  pressure  and  buoyancy  effects  in  early  time  periods  for  air  are  show". 
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MATHEMATICAL  FORMULATION  AND  PHYSICAL  DESCRIPTION 


Consider  the  motion  of  a viscous  fluid  within  a rectangular  enclosure  as 
shown  In  Fig.  1.  The  fluid  Is  Initially  motionless  with  a uniform  temperature 
T0.  The  enclosure  walls  are  also  at  this  temperature.  The  temperature 
difference  Tw-T0  is  initiated  at  time  0 to  Induce  the  flow  within  the 
enclosure.  The  variations  In  all  relevant  physical  properties  are  taken  Into 
account.  However  kinetic  energy.  Internal  heat  sources  and  Irreversible  viscous 
dissipation  In  the  energy  equation  are  not  considered. 

The  governing  equation  for  an  compressible  fluid  with  variable  properties 

will  be  as  follows. (1) 

continuity  equation: 


(1) 


momentum  equations: 


(2) 


where  viscous  forces  on  element  are 
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energy  equation: 
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where 


thermal  expansion  of  the  fluid. 
The  initial  conditions  are: 


is  the  volumetric  coefficient  of 


3P 

3x  - - P9 


(9) 

(10) 

(11) 


and  boundary  conditions  are: 

t > 0 u s v = o at  all  wal Is  (12) 

T * T0  at  all  walls  except  (13) 

heated  surface 

I = Tw  at  heated  surface  (14) 

NUMERICAL  SOLUTION  OF  THE  EQUATION 
Finite  difference  formulation 

The  left  hand  sides  of  eqs.  (1),  (2)  and  (3)  can  be  manipulated  by  noting 

that  for  arbitrary  variable  F as  follows. 
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Using  eq.  (15),  continuity  eq.  (1)  can  be  written 


9(p)  + 9{(p)u)_  + = 0 

3x  9x 


momentum  eq.  (2)  is 

3(pu) ■ 9( (pu)u)  3{(pu)v) 
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and  energy  eq.  (8)  is  also  written  as  follows. 


(19) 
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The  numerical  scheme  shown  In  Fig.  2 employs  a 12x12  rectangular  grid 
system  In  the  x and  y direction  with  a total  of  144  grid  points  and  Ax,  Ay  are 
the  grid  size  in  both  x and  y direction  respectively.  This  numerical  method  is 
based  on  the  simplest  explicit  finite  difference  approximations  to  the  governing 
differential  equations  which  will  be  obtained  at  a finite  number  of  grid  points 
having  coordinate  x s iAx  and  Y * jAY  except  next  to  the  wall  where  the  mesh 
size  is  one-half  the  mesh  within  the  cavity.  All  grid  points  are  evaluated  at 
discrete  times  rn.  The  values  of  all  physical  properties  at  each  grid  point 
should  be  thought  of  as  average  value  over  a small  volume  of  fluid  in  Fig.  2. 

The  finite  difference  approximations  on  the  derivatives  of  the  arbitrary 
variable  F at  the  grid  point  (1,  j)  in  advancing  from  time  rn  to  the  new  level 
Tn+1  = rn  +Ar  may  be  written  as  follows. 


At  (21) 

The  finite  difference  scheme  for  spatial  derivatives  uses  central  differ- 
ences except  next  to  the  wall.  The  first  order  difference  is 
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(22) 


the  second  order  finite  difference  is 
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(23) 


the  cross  finite  difference  is 
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(24) 


d2  F 
dxdy 


(Fi4i,J+i  - Fi-i  .j-n 


- Fn  + Fn 
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4AxAy 


The  nonlinear  convection  terms  cause  the  main  difficulties  to  achieve  a stable 
numerical  method.  For  some  difference  methods,  the  rate  of  heat  removal  may 
differ  from  the  rate  of  heat  addition  at  steady  state.  Torrance  (5)  tested 
several  methods  for  differencing  the  convection  term.  In  this  calculation, 
Torrance's  V method  (5)  is  employed  as  follows. 
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where 
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When  these  approximations  are  Introduced  into  eqs.  (16),  (17),  (18)  and  (19)  we 
obtain 
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where  for  an  ideal  gas 
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Although  the  finite  difference  of  convection  terms  in  eqs.  (29),  (30)  and  (31) 
are  shown  for  only  one  condition  of  eq.  (25),  appropriate  finite  difference  of 
convection  terms  should  be  used  In  actual  computation  work. 

NUMERICAL  PROCEDURE 


The  calculation  proceeds  by  explicitly  advancing  p,  n,  v and  T with 
difference  forms  of  eqs.  (28),  (29),  (30)  and  (31).  Also  pressure  P is 
calculated  explicitly  from  an  equation  of  state  using  P and  T.  Fluid  within 

enclosure  is  initially  at  a uniform  temperature  To  and  at  rest.  Here,  for 
preliminary  studies,  we  consider  a rectangular  enclosure  of  height  Xmax 
(0.1m),  width  Ymax  (0.05m)  and  vertical  heated  wall  (0.04m)  which  is 
located  in  the  middle  part  of  left  side  wall,  as  shown  In  figure  1. 

During  any  one  time-step,  all  values  appearing  in  the  right  side  of  eqs. 

(28),  (29),  (30)  and  (31)  are  treated  as  constants.  In  the  first  place,  the 

new  (pU)?+l  and  (pV)?+!  at  all  Interior  grid  points  may  be  obtained  from 
1 » J • » J n+1 

successive  momentum  eqs.  (29)  and  (30).  Then  new  density  should  be 

calculated  from  continuity  eq.  (28)  substituting  (pU)J+l  and  (pV)"+1.  Into  eq. 
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n+i 

(28).  A new  temperature  Tv  j is  obtained  from  energy  eq.  (31)  using  the  value 

of  just  computed.  Finally  new  velocities  and  V?+l  are  calculated 

• »J  i »j  i 

mathematically  from  (pU)^j  and  (pV)!Jj  using  the  value  of  pn+]  as  follows: 


Un+1 


. (pU) 


n+1 

hi 


>n+! 

i.J 


(32) 


New  pressure  P.  • are  calculated  from  the  equation  of  state  using  new 

1 * J 

temperatures  and  densities  which  are  already  computed  at  Interior  grid  points, 
though  not  at  the  wall.  Pressures  at  the  wall  are  obtained  by  quadratic 
extrapolations.  This  process  is  repeated  in  time,  provided  the  time-step  is 
sufficiently  small.  The  time-step  Axf  has  been  restricted  to  10-5$  or  less 
[5]  in  this  computational  work.  This  value  corresponds  to  less  than  the  time 
interval  for  a sound  wave  to  propagate  across  the  mesh  size  y as  follows 


(33) 


DISCUSSION  OF  RESULTS 

Numerical  calculations  have  been  carried  out  for  air  within  rectangular 
enclosure  (0.1  x 0.05  m2)  in  a time  periods  from  0 to  40  [ms].  The  fluid 
conditions  and  the  imposed  temperature  differences  correspond  to  a Grashof 
number  of  approximately  8.8  x 10^.  Figure  3-(l)  ~ (10)  shows  the  velocity 
vectors  at  Intervals  of  0.03  ms  from  0.03  to  0.30  ms,  where  the  dominant  motion 
Is  normal  to  the  heated  wall  (and  gravity).  The  absence  of  a vector  at  a given 
grid  point  means  that  the  magnitude  of  the  calculated  velocity  at  that  point  was 
less  than  5%  of  the  maximum  velocity  at  any  of  the  grid  points  at  that  instant 
of  time.  The  disturbance-front  separating  the  region  of  non-zero  fluid 
velocities  from  the  region  of  essentially  static  fluid  is  seen  to  move  away  from 
the  heated  wall  with  the  velocity  of  sound.  At  t = 0.15  ms  (figure  3-(5) ),  the 
disturbance  front  reaches  the  right  side  wall.  In  figures  3- (6 ) to  (10),  the 
fluid  motion  becomes  complicated  by  the  sum  of  many  phases  and  amplitudes  of 
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motion  with  multiple  reflection  from  all  walls.  No  gravitational  effects  can  be 
seen. 

Figure  4 — (1)  to  (10)  shows  the  velocity  vectors  at  later  times,  from 
20.03  to  20.30  ms.  The  gravitational  contributions  to  the  fluid  motion,  which 
causes  assymmetry  around  the  horizontal  center  line,  is  still  quite  small 
compared  to  the  motions  induced  by  the  initial  expansion  wave  away  from  the 
heated  surface.  In  order  to  distinguish  the  growth  of  natural  convection,  the 
velocity  components  Ug^2  and  Ug  2 are  shown  in  Fig.  11  for  increasing  time. 

Ug  2and  ug  2are  vertlca^  components  of  velocity  at  grid  points  (5,2)  and 
(8,2)  respectively.  These  grids  points  are  symmetrical  located  about  the 
horizontal  center  line  of  the  heated  surface.  The  data  are  plotted  at  every 
hundreth  time  step,  which  causes  the  apparent  sawtooth  character.  Nevertheless, 
the  superposition  of  various  amplitudes  and  phases,  mentioned  above,  is  clearly 
evident.  The  dotted  lines  Indicate  the  values  of-Ugj2,  i.e.  symmetrical  values 
of  Ug  2about  bhe  zero  velocity.  In  the  early  time  periods,  less  than  about  5 
[ms],  the  velocity  components  Ug  ^ and -Ug s2are  eQual  within  1%  or  better. 

After  that  time  they  begin  to  deviate,  and  this  is  a manifestation  of  buoyancy 
force.  The  difference  between  Ug  2 and  ~U5,2»  illustrated  with  shadow, 
represents  the  growth  of  natural  convection.  The  values  of  the  shadow  of 
deduced  from  Fig.  11  are  shown  in  Fig.  12.  Also  the  difference  of  horizontal 
components  at  both  points  are  shown  in  same  figure  with  the  dotted  curve.  The 
magnitudes  of  velocity  component  at  upper  point  of  (8,2)  are  larger  than  one  at 
lower  point  of  (5,2)  for  either  horizontal  or  vertical  component,  and  this  is 
what  should  be  expected  on  physical  grounds  In  a natural  convection  heat 
transfer.  It  is  of  interest  to  note  that  the  natural  convection  flow  near  the 
heated  wall  induced  by  the  buoyancy  force  develop  continuously  and  smoothly  as 
shown  in  Fig.  12  and  stream  lines,  which  always  close  for  incompressible  fluid 
'flow,  will  not  do  so  in  this  case. 

On  the  other  hand,  the  temperature  field  in  the  vicinity  of  the  heated  wall 
is  essentially  that  of  pure  conduction  for  this  range,  and  isotherms  are 
practically  symmetrical  to  the  heated  wall. 

The  pressure  and  velocity  fluctuation  at  a near-mid  point  of  the  enclosure 
are  shown  in  Figs.  5 and  6 for  time  to  1.0  millisecond.  The  relation  between 
pressure  and  velocity  fluctuation  Is  not  distinct  in  Fig.  5 and  6 for  this 
physical  model  which  has  the  heated  wall  at  the  middle  part  of  the  left  side 
wall  as  in  Fig.  2.  In  order  to  reduce  the  influence  of  reflection  at  upper  and 
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lower  wall,  the  calculation  was  repeated  with  the  heating  surface  extending  over 
the  complete  left  wall  of  the  cavity.  The  calculated  pressure  and  horizontal 

component  of  velocity  are  shown  in  Fig.  7 and  8.  It  is  quite  evident  that  the 
frequency  of  pressure  fluctuation  is  two  times  that  of  the  velocity.  The 
average  pressure  in  the  enclosure  increases  as  heat  is  added,  and  is 
proportional  to  the  heated  area,  as  seen  in  Fig.  5 and  Fig.  7: 


The  results  of  pressure  and  velocity  fluctuation  in  later  time  periods  from 
30.5  to  32.0  [ms]  are  shown  in  Figs.  9 and  10  respectively.  From  these  figures 
we  are  not  directly  able  to  make  clear  the  correlation  between  pressure  and 
velocity  fluctuation.  Therefore  further  study,  using  spectrum  analysis  (or 
something  similar)  should  be  required  to  fully  understand  this  phenomena  where 
the  fluid  motion  is  the  sum  of  many  phases  and  amplitudes  of  motion  with 
multiple  reflection  from  the  walls. 

Numerical  procedures  for  solutions  of  heat  transfer  equations  In  the 
time-dependent  domain  may  fall  Into  two  categories,  explicit  and  Implicit. 

These  two  types  of  difference  equations  have  previously  been  studied  in  which 
explicit  difference  equations  are  simple  to  solve,  but  which  require  a large 

number  of  time  steps  of  limited  size,  and  implicit  difference  equations  do  not 
limit  the  time  step  but  which  do  require  iteration  at  each  time  step  in  the  sol- 
ution. Therefore,  explicit  procedures  are  convenient  under  conditions  where  a 
sufficiently  large  time  step,  consistent  with  computational  stability,  can  be 
used.  In  order  to  examine  the  accuracy  of  this  computational  results,  the 
smaller  time  step  Ar  of  10-4  [ms]  which  is  one  tenth  of  normal  time  step 
have  been  used  for  early  time  periods  from  0 to  0.3  [ms].  The  agreement  with 
the  solutions  of  vy^g  for  this  calculation  of  small  time  step  and  the 
prior  one  Is  better  than  0.1%. 

The  numerical  calculations  have  been  performed  on  a large  digital  computer. 
The  execution  time  at  each  time  step  is  approximately  33  CP  millisecond. 
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Figure  1.  Physical  model  and  coordinates. 
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Figure  2.  Schematic  diagram  of  the  numerical  method. 
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Figure  3-(7).  Velocity  field  at  t-0.21  ms.  Vmax«0.25  mm/s. 
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Figure  3— ( 6 ) . Velocity  field  at  t*0.18  ms.  • 427  mo  s . 
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Figure  3-(5).  Velocity  field  at  t«0.15  ms.  Vmax-0.428  mm/s. 
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Figure  3-(4).  Velocity  field  at  t**0.12  ms.  vmax“0*376  mm/s. 
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Figure  3-(3).  Velocity  field  at  t-0.09  ms.  Vmax-0.566  mm/s. 
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Figure  3— ( 2 ) - Velocity  field  at  t-0.06  ms.  Vmax-0.368  mm/ s . 


21 


(UU)  X 


Figure  3—  Cl).  Velocity  field  at  t-0.03  ms.  vmax"0*521  mm/s. 
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Figure  3-(8).  Velocity  field  at  t«0.24  ms.  Vmax-0.498  mm/s. 
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Figure  3-(9) . Velocity  field  at  t-0.27  ms.  V^x-0.214  mm/s. 
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Figure  3-(10).  Velocity  field  at  t-0.30  ms.  Vmax-0.544  mo's. 
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Figure  4-(l).  Velocity  field  at  t-20.03  ms.  Vmax-0.343  mm/s. 
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Figure  4-(2).  Velocity  field  at  t-20.06  ms.  Vmax**0.380  s . 
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Figure  4-(3).  Velocity  field  at  t-20.09  ms.  mm/s. 
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Figure  4-(4).  Velocity  field  at  t-20.12  ms.  Vmax-0.567  mm's. 
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Figure  4-(5).  Velocity  field  at  t-20.15  ms.  Vn^x-0.339  mm/s. 
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Figure  4-(6).  Velocity  field  at  t-20.18  ms.  Vmax-0.588  mo/s. 
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Figure  4-(7).  Velocity  field  at  t-20.21  ms.  Vmax“0*2^0  mm/8* 
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Figure  4-( 8 ) . Velocity  field  at  t-20.24  ms.  Vmax-0.527  mm  s 
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Figure  4-(9).  Velocity  field  at  t-20.27  ms.  V^-0.325  mm/s. 
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Figure  4-(10).  Velocity  field  at  t-20.30  ms.  Vmax-0.531  mm.'s. 
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Appendix  B 

List  of  variables  [Main  Program] 

These  are  founded  on  computer  program  list  of  1. 


Program  Symbol 

Definition 

A 

=g/(RT),  to  use  initial  density  distribution 

ADATE 

Date  of  computation 

AFPX 

Dummy  argument,  see  sub.  AVEVELO 

AFPY 

t.  .• 

ATIME 

Time  of  computation 

AUPX 

Average  velocity  u(i)  = 0.5[u(i+l,j)  + u(i,j)] 

DUPY 

Average  velocity  u(j)  = 0.5[u(i,j+l)  + u(i,j)] 

AVPX 

Average  velocity  v(i)  * 0.5[v(i+l,j)  + v(i,j)] 

AVPY 

Average  velocity  v(j)  = 0.5[v(i,j+l)  +v(i,j)] 

CODX 

dk/dx 

CODY 

dk/dy 

CONE 

Thermal  conductivity,  k 

CP 

Specific  heat,  Cp 

CPDX 

dCp/dx 

CPDY 

dCp/dy 

DFDX 

Dummy  argument  of  d/dx 

DFDXX 

d2/dx  2 

DFDXY 

" d2/dxdy 

DFDY 

" d/dy 

DFDYY 

" d2/dy2 

DPDR 

dP/dP 

DPUDX 

* d(Pu)/dx  in  energy  equation 

DPVDY 

- d(Pv)/dx 

DRDT 

- dP/dT 

DTDX 

- dT/dx 

DTDXX 

- d2T/dx2 

DTDY 

- dT/dy 

DTDYY 

- d2T/dy2 

DTIME 

-=  At  time  step 
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Line  No 

. Program  Symbol 

Definition 

197 

| DTUDX 

8(Tu)/8x  difference  of  convection  term 

197 

| DTVDY 

- 

8(Tv)/dx 

32 

| DTW 

s 

Tw-T0,  temperature  difference 

| DUDX 

c 

du/ dx 

| DUDXX 

m 

82u/dx2 

| DUDXY 

sr 

82u/dx8y 

| DUDY 

s 

8u/dy 

| DUDYY 

ss 

82u/8y2 

| DURDX 

= 

8(Pu)/dx 

| DURDY 

ac 

8(Pu)/8y 

| DVDX 

= 

8v/dx 

| DVDXX 

s 

d2v/dx2 

| DVDXY 

s 

82v/ dxdy 

| DVDY 

ss 

8v/8y 

| DVDYY 

= 

82v/8y2 

| DVRDX 

* 

8 (Pv)/dx 

| DVRDX 

s 

8( Pv)/8y 

1 DX 

- 

dx 

| DY 

- 

dy 

1 G 

- 

g,  acceleration  of  gravity 

| I 

i 

th  grid  point  in  x direction 

| ICO 

- 

0,  printout  at  first  time  step 

74,225 

| ICOUNT 

counter  on  the  number  of  time-step  for  printout 

232 

1 « 

integer  of  changing  of  I th  order  for  printout 

59 

1 IS1 

The  lowest  grid  point  of  heated  surface 

59 

I IS2 

The  highest  grid  point  of  heated  surface 

225 

| I WRITE 

Controller  integer  value  of  printout 

1 J 

j 

th  grid  point  in  y direction 

1 K 

K 

has  1,2  and  3.  K=3 : newest  one;  K=1 : old  one 

1 M 

M 

th  grid  point  is  correspond  to  Xmav 

| Ml 

- 

M-l 

| M2 

ST 

M-2 

1 N 

N 

th  grid  point  corresponds  to  ymax 

N1 

- 

N-l 

| N2 

m 

N-2 
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Line  No. 

Program  Symbol 

Definition 

PRDX 

- 

PRDY 

- 

PRESS 

pressure,  P 

R 

density,  P 

149 

RU 

* Pu 

RUDX 

■ d(Pu)/dx 

157 

RV 

“ PV 

RVDY 

■ d(pv)/d  y 

T 

temperature 

32 

TB 

static  temperature  or  wall  temp,  except  heated  wall 

269 

TMAX 

limitation  of  computing  time 

61 

TW 

* Tw,  temperature  at  heated  wall 

73 

TYME 

* t + At  , increment  time  in  computation 

199 

T1 

some  term  in  energy  equation 

200 

T2 

• 0 

201 

T3 

•« 

202 

T4 

•t 

204 

T5 

•• 

205 

T6 

M 

206 

T7 

M 

U 

velocity  component  in  x direction 

UDV 

is  not  used,  only  dimension 

UVMG 

is  not  used,  only  dimension 

143 

U1 

some  term  in  mementum  equation 

144 

U2 

•t 

145 

U3 

M 

147 

U4 

” 

V 

velocity  component  in  y direction 

145 

VIDX 

= dp/dx 

147 

VIDY 

= d^/dy 

vise 

* P,  viscosity 

151 

VI 

some  term  in  momentum  equation 

152 

V2 

" 

153 

V3 

M 

155 

V4 

M 

45 

X 

* X 

30 

XMAX 

“ xmax 

31 

YMAX 

■ ymax 
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[Subroutine] 


Line 

*** 


*** 


*** 


*** 


*** 


*** 


No. 

Program  Symbol 

Definition 

PROPER 

Subprogram  Name  [calculation  for  properties] 

COND 

thermal  conductivity  " k 

CP 

specific  heat,  « Cp 

DPDR 

- dP/dP 

DPDT 

- dP/dT 

DRDT 

- dP/d T 

M3 

- M-3 

N3 

- N-3 

PRESS 

pressure,  ■ P 

R 

density,  « P 

T 

temperature,  * T 

vise 

absolute  viscosity,  - m 

FIRSTDF 

Subprogram  Name  [difference  for  first  order] 

DFDX 

F is  dummy  value,  * dF/dx 

DFDY 

“ dF/dy 

DX 

mesh  size  in  x direction 

DY 

mesh  size  in  y direction 

F 

dummy  argument  value 

SECONDDF 

DFDXX 

Subprogram  Name  [difference  for  second  order] 
- d^F/tfx2 

DFDYY 

- d2F/ (9y 2 

DXX 

- (Ax)2 

DYY 

- (AY)2 

QDXX 

* 1/4  (Ax)2 

QDYY 

- 1/4  (AY)2 

CROSSDF 

Subprogram  Name  [cross  difference] 

DXY 

* Ax  * Ay 

CONVDF 

Subprogram  Name  [difference  for  convection  term] 

AUPX 

see  page  1 of  this  appendix 

AUPY 

AVPX 

AVPY 

DUFDX 

* d(uF)/dx  F is  dummy  argument 

DVFDY 

■ d(vF)/dy 

AVEVELO 

Subprogram  Name  [calculation  of  average  velocity] 

AFPX 

- AUPX  or  AVPX 

AFPY 

• AUPY  or  AVPY 
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